// Matrix.h: interface for the CMatrix class.
// Part of KenTools library.
// Written by Kwan-Yee Kenneth Wong 2000.
//
//////////////////////////////////////////////////////////////////////
#pragma once
// define MATRIX_THROW_EXCEPTIONS to have exceptions thrown on invalid requests
// defined by default for debug mode
#ifdef _DEBUG
#   define MATRIX_THROW_EXCEPTIONS
#endif

#ifdef MATRIX_THROW_EXCEPTIONS
/** if you have boost available, define MATRIX_USE_BOOST_EXCEPTIONS for
    more detailed error info. If you do not want to depend on Boost, simply
    don't define, all the additional error info is then stripped by the
    PreProcessor so we don't get any compilation errors.
    Catch boost exceptions with:
    catch (boost::exception& e)
    {
        cout << boost::diagnostic_information(e) << endl;
    }
    as a boost exception derives from std::exception, it can also be caught
    with the below, at the cost of loosing the extra error information.
    catch (exception& e)
    {
        cout << e.what() << endl;
    }
*/
#   define MATRIX_USE_BOOST_EXCEPTIONS
#   include <exception>
#   include <string>
#   ifdef MATRIX_USE_BOOST_EXCEPTIONS
#       pragma warning(push)
#       pragma warning(disable: 4996)
#       include <boost/exception/all.hpp>
#       include <boost/tuple/tuple.hpp>
#       define MATRIX_THROW_EXCEPT1(e) \
            BOOST_THROW_EXCEPTION(boost::enable_error_info(std::exception(e)))
#       define MATRIX_THROW_EXCEPT2(e,i) \
            BOOST_THROW_EXCEPTION(boost::enable_error_info(std::exception(e)) << i)

        typedef boost::error_info<struct excep_max_idx_     ,int>           excep_max_idx;
        typedef boost::error_info<struct excep_this_index_  ,int>           excep_this_index;
        typedef boost::error_info<struct excep_attrib_      ,std::string>   excep_attrib;

        // matrix sizes
        typedef boost::error_info<struct excep_this_nrows   ,int>           excep_nrows_t;
        typedef boost::error_info<struct excep_this_ncols   ,int>           excep_ncols_t;
        typedef boost::tuple     <excep_nrows_t,excep_ncols_t>              excep_size_this;
        typedef boost::error_info<struct excep_other_nrows  ,int>           excep_nrows_o;
        typedef boost::error_info<struct excep_other_ncols  ,int>           excep_ncols_o;
        typedef boost::tuple     <excep_nrows_o,excep_ncols_o>              excep_size_other;
#       pragma warning(pop)
#   else
#       define MATRIX_THROW_EXCEPT1(e) \
            throw std::exception(e)
#       define MATRIX_THROW_EXCEPT2(e,i) \
            throw std::exception(e)
#   endif
#endif

#include <vector>

template <class Type> 
class CMatrix
{

public:
    // constructors
    CMatrix();
    CMatrix(int r, int c);
    CMatrix(int r, int c, Type v);
    CMatrix(const CMatrix<Type>& V);        // from vector
    CMatrix(const std::vector<Type>& M);    // from matrix
    CMatrix(int Len);	// as vector

    // destructors
    virtual ~CMatrix();

    // initializations
        // if the size of the current allocated memory is the same as requested,
        // SetSize will only change the values of m_bRows and m_nCols and all the
        // values in the original matrix are retained

        // set the size of the matrix
    virtual bool SetSize(int Rows, int Cols);
        // set the size of the matrix and initialize every element to Val
    virtual bool SetSize(int Rows, int Cols, Type Val);
        // set the size of the vector
    virtual bool SetSize(int Len);
        // set every element to zero
    virtual void SetZeros();

    // retrievers
        // return number of rows
    inline int Rows() const {return m_nRows;}
        // return number of columns
    inline int Cols() const {return m_nCols;}
        // return the length of the vector
    inline int Len() const {return m_nRows*m_nCols;}
        // check if the matrix is valid
    inline bool isValid() const {return m_nRows!=0;}
        // return a pointer to a row
    inline const Type* Row(int r) const {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (r>=Rows())
            MATRIX_THROW_EXCEPT2("Matrix const Row: Provided bad index into matrix",
                excep_size_this(Rows(),Cols()) <<
                excep_max_idx(Rows()-1) <<
                excep_this_index(r));
#endif
        return (m_pData + r*m_nCols);
    }
    inline Type* Row(int r) {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (r>=Rows())
            MATRIX_THROW_EXCEPT2("Matrix Row: Provided bad index into matrix",
                excep_size_this(Rows(),Cols()) <<
                excep_max_idx(Rows()-1) <<
                excep_this_index(r));
#endif
        return (m_pData + r*m_nCols);
    }
        // return a pointer to the internal buffer
    inline const Type* Buffer() const {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (!isValid())
            MATRIX_THROW_EXCEPT1("Matrix const Buffer: Empty matrix");
#endif
        return m_pData;
    }
    inline Type* Buffer() {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (!isValid())
            MATRIX_THROW_EXCEPT1("Matrix Buffer: Empty matrix");
#endif
        return m_pData;
    }
        // retrieve the sub-matrix
    virtual bool SubMat(int r, int rr, int c, int cc, CMatrix<Type>& M) const;
    virtual CMatrix<Type> SubMat(int r, int rr, int c, int cc) const;
        // return a row as a separate vector
    virtual bool GetRow(int r, CMatrix<Type>& M) const;
    virtual CMatrix<Type> GetRow(int r) const;
        // return a column as a separate vector
    virtual bool GetCol(int c, CMatrix<Type>& M) const;
    virtual CMatrix<Type> GetCol(int c) const;
        // return the vector norm
    virtual double GetVecNorm() const;
        // return the homogeneous vector
    virtual CMatrix<Type> Homogeneous(Type lastdim = 1);
    virtual void Normalize();

    // manipulators
        // check if the index is okay
    virtual bool isInside(int r, int c) const;
        // return the trace (sum of the diagonal elements)
    virtual Type Trace() const;
        // perform QR decomposition
    virtual bool QR(CMatrix<Type>& Q, CMatrix<Type>& R) const;
        // perform modified QR decomposition such that R is lower-triangular
    virtual bool QR2(CMatrix<Type>& Q, CMatrix<Type>& R) const;
        // perform Cholesky decomposition where L is lower-triangular and M = L.L^T 
        // only the upper triangular part of M will be considered and used in the decomposition
    virtual bool Chol(CMatrix<Type>& L) const;
        // perform singular value decomposition
    virtual bool SVD(CMatrix<Type>& U, CMatrix<Type>& D, CMatrix<Type>& V) const;
        // perform singular value decomposition with descending singular values in D
    virtual bool SVD2(CMatrix<Type>& U, CMatrix<Type>& D, CMatrix<Type>& V) const;
        // compute rank-n approximation with singular value matrix.
    virtual void GetRankNApproximation(const int& nRank_, CMatrix<Type>* pSingularValues_,CMatrix<Type>* pRank3Approximation_) const;
        // calculate the Frobenius Norm with 
    virtual double GetFrobeniusNorm(CMatrix<Type>& M_) const;
        // return the determinant of the matrix
    virtual double Det() const;
        // return the transpose
    virtual bool Transpose(CMatrix<Type>& M) const;
    virtual CMatrix<Type> Transpose() const;
        // return the inverse
    virtual bool Inverse(CMatrix<Type>& M) const;
    virtual CMatrix<Type> Inverse() const;
        // solve for X where M*X = B
    virtual bool Solve(const CMatrix<Type>& B, CMatrix<Type>& X) const;
    virtual CMatrix<Type> Solve(const CMatrix<Type>& B) const;
        // form a diagonal matrix assuming the values from matrix M row-wise,
        // failure will clean up this matrix
    virtual CMatrix<Type>& Diag(const CMatrix<Type>& M);
        // form a diagonal matrix assuming the values from array M with length l,
        // failure will clean up this matrix
    virtual CMatrix<Type>& Diag(const Type M[], int l);
        // form a diagonal matrix assuming the same value V,
        // the size of the matrix must be set beforehand
    virtual CMatrix<Type>& Diag(Type V);
        // form a vector from the diagonal entries of the matrix
    virtual CMatrix<Type> Diag() const;
        // set values of the submatrix to M,
        // failure will keep this matrix untouched
    virtual bool Set(int r, int c, const CMatrix<Type>& M);
        // set values of the submatrix to M(r1:r2,c1:c2)
        // failure will keep this matrix untouched
    virtual bool Set(int r, int c, const CMatrix<Type>& M, int r1, int r2, int c1, int c2);
        // set values of the submatrix to T,
        // failure will keep this matrix untouched
    virtual bool Set(int r, int rr, int c, int cc, const Type& T);
        // set the row of M,
        // failure will keep this matrix untouched
    virtual bool SetRow(int r, const CMatrix<Type>& Row);
        // set the column of M
        // failure will keep this matrix untouched
    virtual bool SetCol(int c, const CMatrix<Type>& Col);
        // swap rows r and rr
        // failure will keep this matrix untouched
    virtual bool SwapRows(int r, int rr);
        // swap columns c and cc
        // failure will keep this matrix untouched
    virtual bool SwapCols(int c, int cc);
        // swap values i and j
        // failure will keep this matrix untouched
    virtual bool Swap(int i, int j);
        // repeat matrix by rxc
    virtual bool RepMat(int r, int c, CMatrix<Type>& M) const;
    virtual CMatrix<Type> RepMat(int r, int c) const;
        // for vector cross product: M.Tx()*V = M x V, with both M and V are 3x1 vectors
    virtual bool Tx(CMatrix<Type>& M) const;
    virtual CMatrix<Type> Tx() const;
        // convolution and polynomial multiplication, treat all matrices as vectors
    virtual bool Conv(const CMatrix<Type>& M, CMatrix<Type>& Result) const;
    virtual CMatrix<Type> Conv(const CMatrix<Type>& M) const;

        // normalized cross correlation, get aligned with the mean
    virtual Type NCC(const CMatrix<Type>& M_) const;
        // normalized cross correlation, get aligned with 0; 
    virtual Type NCC2(const CMatrix<Type>& M_) const;
    // array multiplication (element-wise)
    virtual bool Times(const CMatrix<Type>& M, CMatrix<Type>& Result) const;
    virtual CMatrix<Type> Times(const CMatrix<Type>& M) const;
        // array division (element-wise)
    virtual bool Div(const CMatrix<Type>& M, CMatrix<Type>& Result) const;
    virtual CMatrix<Type> Div(const CMatrix<Type>& M) const;
        // calculate the mean for each row
    virtual bool RowMean(CMatrix<Type>& M) const;
    virtual CMatrix<Type> RowMean() const;
        // calculate the mean for each column
    virtual bool ColMean(CMatrix<Type>& M) const;
    virtual CMatrix<Type> ColMean() const;
        // calculate the mean for the whole matrix
    virtual double Mean() const;
    virtual double Max() const;
    virtual double Max(int *pc,int *pr) const;
    virtual double Min() const;
        // clean up the memory and set size to zero
    virtual void CleanUp();

    // operators
        // access M(r,c) without range checking
    inline const Type& operator () (int r, int c) const {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (r>=Rows() || c>=Cols())
            MATRIX_THROW_EXCEPT2("Matrix const M(r,c) operator: Provided bad indices into matrix",
                excep_size_this(Rows(),Cols()) <<
                excep_attrib(r>=Rows()? "row" : "col") <<
                excep_max_idx(r>=Rows()? Rows()-1 : Cols()-1) <<
                excep_this_index(r>=Rows()? r : c));
#endif
        return m_pData[r*m_nCols+c];
    }
    inline Type& operator () (int r, int c) {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (r>=Rows() || c>=Cols())
            MATRIX_THROW_EXCEPT2("Matrix M(r,c) operator: Provided bad indices into matrix",
                excep_size_this(Rows(),Cols()) <<
                excep_attrib(r>=Rows()? "row" : "col") <<
                excep_max_idx(r>=Rows()? Rows()-1 : Cols()-1) <<
                excep_this_index(r>=Rows()? r : c));
#endif
        return m_pData[r*m_nCols+c];
    }
        // access M(i) as a vector without range checking
    inline const Type& operator () (int i) const {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (i >= Len())
            MATRIX_THROW_EXCEPT2("Matrix const M(i) operator: Provided bad index into matrix",
                excep_size_this(Rows(),Cols()) <<
                excep_max_idx(Len()-1) <<
                excep_this_index(i));
#endif
        return m_pData[i];
    }
    inline Type& operator () (int i) {
#ifdef MATRIX_THROW_EXCEPTIONS
        if (i >= Len())
            MATRIX_THROW_EXCEPT2("Matrix M(i) operator: Provided bad index into matrix",
                excep_size_this(Rows(),Cols()) <<
                excep_max_idx(Len()-1) <<
                excep_this_index(i));
#endif
        return m_pData[i];
    }
        // return the submatrix (r:rr,c:cc)
    virtual CMatrix<Type> operator () (int r, int rr, int c, int cc) const;

    // conversion operator
    //operator             Type () const; // if matrix is scalar, extracts the scalar
    operator std::vector<Type>() const; // if matrix is vector, extracts the vector

        // negation
    virtual CMatrix<Type> operator - () const;

        // arithmetic operations with matrix M
    virtual CMatrix<Type> operator + (const CMatrix<Type>& M) const;
    virtual CMatrix<Type> operator - (const CMatrix<Type>& M) const;
    virtual CMatrix<Type> operator * (const CMatrix<Type>& M) const;
        
        // arithmetic operations with scalar T
    virtual CMatrix<Type> operator + (const Type& T) const;
    virtual CMatrix<Type> operator - (const Type& T) const;
    virtual CMatrix<Type> operator * (const Type& T) const;
    virtual CMatrix<Type> operator / (const Type& T) const;
        
        // assignment operations with matrix arithmetics (all prefix const eliminated by jerry)
    virtual CMatrix<Type>& operator  = (const CMatrix<Type>& M);
    virtual CMatrix<Type>& operator += (const CMatrix<Type>& M);
    virtual CMatrix<Type>& operator -= (const CMatrix<Type>& M);
    virtual CMatrix<Type>& operator *= (const CMatrix<Type>& M);
        
        // assignment operations with scalar arithmetics (all prefix const eliminated by jerry)
    virtual CMatrix<Type>& operator += (const Type& T);
    virtual CMatrix<Type>& operator -= (const Type& T);	
    virtual CMatrix<Type>& operator *= (const Type& T);
    virtual CMatrix<Type>& operator /= (const Type& T);

    virtual CMatrix<Type>& operator = (const std::vector<Type>& V);             // assignment operator

    // serializations
    //virtual void SerializeTxt(CArchive& ar);
    //virtual void Serialize(CArchive& ar);
    //
    //// debug purpose
    //virtual void DumpMatrix();

protected:
    int   m_nRows;
    int   m_nCols;
    Type* m_pData;

    // helper functions
    double Pythag(double a,double b) const;
};





template <class Type>
CMatrix<Type> operator + (const Type& T, const CMatrix<Type>& M);

template <class Type>
CMatrix<Type> operator - (const Type& T, const CMatrix<Type>& M);

template <class Type>
CMatrix<Type> operator * (const Type& T, const CMatrix<Type>& M);

template <class Type>
CMatrix<Type> operator * (const CMatrix<Type>& M, const Type& T);


template <class Type>
CMatrix<Type> operator / (const Type& T, const CMatrix<Type>& M);

template <class Type>
Type dot(const CMatrix<Type>& A, const CMatrix<Type>& B);

template <class Type>
CMatrix<Type> cross(const CMatrix<Type>& A, const CMatrix<Type>& B);

// macro for instantiating different types of the matrix class
#define Instantiate_Matrix(Type)\
    template class CMatrix<Type>;\
    template CMatrix<Type> operator + (const Type& T, const CMatrix<Type>& M);\
    template CMatrix<Type> operator - (const Type& T, const CMatrix<Type>& M);\
    template CMatrix<Type> operator * (const Type& T, const CMatrix<Type>& M);\
    template CMatrix<Type> operator / (const Type& T, const CMatrix<Type>& M);\
    template Type dot(const CMatrix<Type>& A, const CMatrix<Type>& B);\
    template CMatrix<Type> cross(const CMatrix<Type>& A, const CMatrix<Type>& B);

typedef CMatrix<double> CDMatrix;



//// program option parser for CMatrix
//#ifndef BOOST_PROGRAM_OPTIONS_NO_LIB
//    #define BOOST_PROGRAM_OPTIONS_NO_LIB
//#endif
//#include <boost/throw_exception.hpp>
//#include <boost/algorithm/string.hpp>
//#include <boost/any.hpp>
//#include <boost/program_options/errors.hpp>
//
//namespace boost { namespace program_options { 
//
//    template<class T, class charT>
//    void validate(boost::any& v, 
//                  const std::vector<std::basic_string<charT> >& s, 
//                  CMatrix<T>*,
//                  int)
//    {
//        if (v.empty()) {
//            v = boost::any(CMatrix<T>());
//        }
//        CMatrix<T>* tv = boost::any_cast< CMatrix<T> >(&v);
//        assert(NULL != tv);
//        for (unsigned i = 0; i < s.size(); ++i)
//        {
//            try {
//                /* We call validate so that if user provided
//                a validator for class T, we use it even
//                when parsing CMatrix<T>.  */
//
//                vector<std::basic_string<charT>> vsRow;
//                vector<std::basic_string<charT>> vsColumn;
//                int nNC;    // Number of columns for matrix
//
//                // test if vector notation is used in input
//                assert(*s[i].begin() == '(' && *s[i].rbegin() == ')');
//
//                // split rows of matrix
//                split( vsRow, s[i].substr(1, s[i].size()-2), is_any_of(";"), token_compress_on );
//
//                // for each row, get the values and put them in matrix
//                for (unsigned kr = 0; kr < vsRow.size(); kr++)
//                {
//                    vsColumn.clear();
//                    split( vsColumn, trim_copy(vsRow[kr]), is_any_of(" "), token_compress_on );
//
//                    if (kr==0) // matrix has size 0 now, set to expected size
//                    {
//                        tv->SetSize(vsRow.size(),vsColumn.size(),0);
//                        nNC = vsColumn.size();
//                    }
//                    else // check if current row has same number of columns as previous
//                        assert(vsColumn.size()==nNC);
//
//                    // validate option values
//                    for (unsigned kc = 0; kc < vsColumn.size(); kc++)
//                    {
//                        boost::any a;
//                        std::vector<std::basic_string<charT> > v;
//                        v.push_back(vsColumn[kc]);
//                        validate(a, v, (T*)0, 0);                
//                        (*tv)(kr,kc) = boost::any_cast<T>(a);
//                    }
//                }
//            }
//            catch (const bad_lexical_cast&) {
//                boost::throw_exception(invalid_option_value(s[i]));
//            }
//        }
//    }
//
//}}